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Belief propagation — a powerful heuristic method to solve inference problems involving a large 
number of random variables — was recently generalized to quantum theory. Like its classical 
counterpart, this algorithm is exact on trees when the appropriate independence conditions are 
met and is expected to provide reliable approximations when operated on loopy graphs. In this 
paper, we benchmark the performances of loopy quantum belief propagation (QBP) in the context 
of finite-temperature quantum many-body physics. Our results indicate that QBP provides reliable 
estimates of the high-temperature correlation function when the typical loop size in the graph is 
large. As such, it is suitable e.g. for the study of quantum spin glasses on Bethe lattices and the 
decoding of sparse quantum error correction codes. 
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I. INTRODUCTION 

Belief propagation is a powerful algorithm designed to 
solve inference problems involving a large number of ran- 
dom variables. It operates on graphical models, where 
variables are located at the vertices of a graph and edges 
encode dependence relations between the variables. The 
algorithm is exact when the underlying graph is a tree, 
but most importantly, it performs remarkably well in cir- 
cumstances where it is not proven to converge to the 
exact solutions, i.e. when the graphical model contains 
loops. It is also highly parallelizablc in the sense that 
each random variable can be associated with a different 
processor, and messages are exchanged between proces- 
sors that are joined by an edge [U, UH, 

These features have made belief propagation an impor- 
tant tool in numerous scientific and technological fields 
ranging from information theory to image recognition, 
and from artificial intelligence to statistical physics. In- 
deed, it is one of the most powerful heuristic algorithms 
to solve problems such as decodingof low-density and 
turbo error correction codes [1, 0, determining the 
phase diagram of quenched disordered systems [1, [1] , and 
random satisfiability problems 0, [13] ■ 

Recently, belief propagation and graphical models were 
generalized to the quantum setting [ll|, • In this arti- 
cle, we characterize the performances of QBP when used 
as a heuristic algorithm to solve inference problems — 
e.g. compute correlation functions — in the context of 
finite-temperature quantum many-body physics. 



II. GRAPHICAL MODELS 

We consider quantum graphical models (G, p) that 
consist of a graph G and an n-bifactor state p. The 
graph G = (V, E) has a set of vertices V and a set of 
edges E. Each w G is a quantum system, with Hilbert 
space TLtj. A n-bifactor state p is a positive operator on 



Ti = (S'tiey'^tn that can be expressed as 

tiGV {u,v)eE 

where Z is some normalization factor, and p^ and 
are positive operators on 7i„ and 7i„ ® Hv , respectively. 
The operators Vu-.v are required to mutually commute 
when n is finite. The product is defined as X •k^^'> 
Y = [X^Y~ X^]" , and has the property of producing a 
positive operator when both X and Y are positive. This 
product is non-commutative except in the limit n —^ oo, 
which defines the product: XqY = lim„^oo = 
g(iogj>f+iogy)^ Both products and reduce to normal 
matrix product when X and Y commute. 

A generic inference problem on a graphical model is 
to compute the reduced density operator on a subset 
W d V of the quantum systems conditioned on the fact 
that a measurement was performed on a disjoint subset 
U C V, where both W and U are of constant size. This 
problem turns out to be equivalent to the seemingly sim- 
pler problem of computing the reduced state on any sub- 
set W of constant size, i.e. pw = Tt^v-w{p}^ where Trx 
denotes the partial trace over a systems in set X. With- 
out additional assumptions on the structure of p, solving 
this problem requires resources that grow exponentially 
with the number of quantum systems \V\. However, the 
solution can sometimes be obtained or approximated by 
QBP in a time polynomial in \V\. 



III. QBP ALGORITHM 

Given a graphical model (G, p), QBP consists of a se- 
quence of exchanges of operator- valued messages between 
neighboring vertices, which carry information about the 
state at other locations in the graph. More precisely, for 
(u, v) S E, the message passed from vertex u to vertex v 
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at time t is an operator on Tiy given by 
ruu^vit) oc Tr„|At„*(") (i^u-.vQQ m„,^„(<-l))| (2) 

where n{u) denotes the neighbors of u. The proportion- 
ahty factor can be chosen so that Trlm^^^} = 1, and 
the messages are initiahzed to„^i,(0) = /. At time t the 
behef buv (t) — which is meant to represent some approx- 
imation of the state puv = Try_,j„{p} for (u, u) G -E — 
is given by 

(3) 

'W^n{u)~v y^n{v) — u 

where all messages are taken at time t. When all oper- 
ators defining the bifactor state commute, QBP reduces 
to the standard belief propagation algorithm [l|, H, d, 0] • 
Since the message update rule Eq. ^ at vertex u de- 
pends only on the incoming messages at that vertex, the 
algorithm can be operated in a highly parallel fashion 
where each quantum system u is associated with a pro- 
cessor, and messages are exchanged between processors 
u and V iff (u, v) G E. Similarly, the beliefs on the pair 
(w, v), Eq. ([3]), can be computed by combining the mes- 
sages received at those vertices. 



n = 1, basic algebra implies that cx Trt,{(/z„^u,) ★^^^ 
(Tr„{/i„ -k^^^ Vu:v} Vv■.w)\^ the operations Tr„ and 
commute so to say. The computation of can thus be 
broken into two steps: z) Compute to„^i, — Tr„{^t, ik-^^^ 

}; a) Compute p-a, cx Tr„{(/Lt„/Zu,)*(^) ( 
When n ^ cxo on the other hand, the operations Tr„ and 
do not commute in general, but they do precisely when 
I(u : w\v) = [l^]. QBP is based on a generalization of 
these observations to arbitrary graphs. 

QBP does not rely on the vanishing of the normalized 
connected correlation functions C{<7a, o'c) = i'^AO'c) ~ 
{aA){crc) or equivalently on the vanishing of 

the mutual information J(yl : C) = S{A)+S{C)-S{AC). 
In many systems, the mutual information is not a priori 
short range. For instance in the T — > limit, the ther- 
mal state of the 1-d Ising model in zero transverse field 
is an equal mixture of all spins up and all spins down, 
which has I{A : B) ~ 1 between any two disjoint regions, 
whereas I{A : C\B) — for any three disjoint regions. 
To compute thermodynamical quantities, one generally 
introduces a symmetry-breaking field that randomly sin- 
gles out either the all-up or all-down state, which both 
have I{A : B) ~ 0. Symmetry-breaking can be a delicate 
issue — for instance, on Cayley trees where a constant 
fraction of vertices live on the boundary [1] — and is 
circumvented by QBP. 



A. Convergence 

In it was shown that when G is a tree and (G, p) 
is either i) a 1-bifactor state [c.f. Eq. ^ with n = 1] 
or ii) a quantum Markov network, QBP yields the exact 
solution in a time proportional to the graph's diameter 
— i.e. buv(t) = Puv for t > diameter (G). Intuitively, this 
means that the algorithm must run for a time sufficiently 
long to allow messages to travel between any pair of ver- 
tices. When operated on loopy graphs, the beliefs do not 
necessarily converge to the correct density operators. A 
good heuristic in that case is to halt the algorithm when 
buv{t) become almost time-independent, which also hap- 
pens in a time roughly equal to the graph's diameter in 
all the models we have investigated. 

A graphical model (G, p) is a quantum Markov net- 
work when the conditional independence conditions I{U : 
{V - n{U) - U)\n{U)) = are met for all U ^V. The 
quantity /(yl : B\C) = S{AC)+S{BC)-S{C)-S{ABC) 
is the quantum conditional mutual information [3; ) 
and S{A) = Tr{/5^ log^PA} is the von Neumann entropy. 
As explained in p^. Il5j|. the vanishing of I {A : C\B) 
is equivalent to the condition pabc — P'b © Pab 
Pbc- This equality is not verified in general, and the 
KuUback-Leibler distance between the right- and left- 
hand side is precisely the conditional mutual information 
T^{pABc\\p~B Pab(z) pbc) = ■■ C\B). 

To understand the workings of QBP, consider a bifac- 
tor state Puvw on the line u — v — w. The reduced state 
on ui is pu, cx TiuviiP'uPvPw) {i^u-.v i^v.w)}- When 



IV. QBP FOR QUANTUM MANY-BODY 

In the context of quantum many-body physics, the in- 
ference problem consists of computing correlation func- 
tions for the thermal state of a system of interacting par- 
ticles. Given a graph G — {V, E), we consider the generic 
Hamiltonian 

veV {u,v)eE 

The thermal state at inverse temperature f3 = 1/T is 
given by p = -^e~^^ where Z — Tr{e~''^} is the par- 
tition function. Defining /x„ = ^-0^^ and Vu-.v = e~^^"^ 
enables us to express any such thermal state as an oo- 
bifactor state, cf. Eq. dT]). 

Despite the fact that thermal states are bifactor states, 
the result from [13] cited above does not imply that cor- 
relation functions can be evaluated exactly and efficiently 
with QBP. This is primarily because G is not necessar- 
ily a tree, but also because thermal states are neither 1- 
bifactor nor quantum Markov networks in general. There 
is no general remedy to the first hurdle, unless the loops 
happen to be very small and can be eliminated by merg- 
ing some vertices. Thus, QBP will need to be executed 
on a loopy graph and it is the primary goal of this pa- 
per to determine the effects of such loops on the per- 
formance of QBP. Two pragmatic solutions, named the 
replica method and sliding window QBP, have been pro- 
posed to overcome the second set of obstacles jl^]. 
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A. Replica 

The general idea of the rephca method is to approx- 
imate the thermal state by a 1-bifactor state on which 
QBP can be executed directly and is guaranteed to con- 
verge in the absence of loops. In a first step, a Trotter- 
Suzuki (TS) decomposition is used to approximate a ther- 
mal state by an 7VT--bifactor state with finite Nr- This 
produces a systematic error that scales as P/N-r- Then, 
in a fashion reminiscent of the replica trick used in the 
study of spin glass, the A'r-bifactor state is replaced by a 
1-bifactor state at the expense of substituting the quan- 
tum system at each vertex by replicas: 



Ml) 



and 



(5) 



where T^^^' is the operator that cyclicly permutes the 
Nt replicas of v. The operators v^-v — e^^^^'-' do not 
commute in general, but this can be fixed in practice on 
sparse graphs by merging some vertices. On a tree, the 
TS decomposition is the only source of error, so accuracy 
e can be achieved at a computational cost that is expo- 
nential in P/e. This method is particularly useful as it 
allows for a direct computation of correlation functions 
at arbitrary distances, see [T^ . 



B. Sliding window 

While all quantum Markov networks are thermal states 
of some local Hamiltonian on G [12] , the converse is not 
true in general. Sliding window QBP is motivated by 
the fact that quantum Markov networks are fixed points 
of coarse graining procedures. Thermal states, regarded 
as cxD-bifactor states, are used directly to implement the 
message passing rule in Eq. ^ with n — oo, except 
that messages are computed not just using the nearest 
neighbors but with all vertices within a distance < £. 
On a line, for instance, vertex j receives a message from 
(j — 1, j — 2, . . .j — t) and one from (j + I,j + 2, . . . ,j + £). 
In that case, sliding window QBP produces the exact so- 
lution efhciently if the conditional mutual information 
dies off at a finite distance. 



V. NUMERICAL RESULTS 

We have numerically implemented the QBP algorithm 
on various graphs for the Ising and Heisenberg model 
whose Hamiltonians are 



Hi ^^g-a^ 

vev 



JuvKcTy and Hh = g„ ■ 



{u,v)GE 



(6) 

respectively, and a — {a^ ,a^) are the usual Pauh 
matrices normalized so that — 11/4. On a line, the 
homogeneous {Juv = 1) Ising model has a zero tem- 
perature phase transition at the critical transverse field 




TEBD 



Sliding window 
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FIG. 1: (Color online) Critical Ising model on infinite line. 
Energy density estimate E using the method of replicas with 
= 10, sliding window QBP with £ = 6, and TEBD with 
X = 150, compared to the exact energy density E obtained 
from fermionization. 



g = (^,0,0). Most of our simulations were performed 
at this critical value, as it is expected to represent the 
"hardest case". Unless otherwise specified, it is hence- 
forth assumed that g ~ (i, 0, 0) and Juv — 1- 

We used QBP to compute the energy density of the 
Ising model on an infinite line. This model can be solved 
exactly by means of a Jordan- Wigner transform that 
maps the interacting spin chain to a collection of free 
fermions [l^ . Figure [1] shows the difference between the 
energy density computed with QBP and its exact value as 
a function of inverse temperature. Also shown are results 
obtained from a superoperator version of time-evolving 
block decimation (TEBD), which combines ideas from 
(20I [2H . Since a line is a tree, the error in the results 
obtained from the replica method is entirely caused by 
the TS decomposition. The results obtained for sliding 
window QBP are in remarkably good agreement with the 
exact value, and can be systematically improved by in- 
creasing £. This reflects the fact that correlations are 
short-ranged in finite-temperature ID models. As ex- 
pected the agreement improves for non-critical g. Re- 
sults obtained for the Heisenberg model on the infinite 
line (not shown) are similar in all aspects. 

To characterize the performance of QBP on more gen- 
eral graphs, we restrict our attention to systems with less 
than 12 spins, allowing comparison to direct brute-force 
numerical solutions. Figure [2] shows the correlation func- 
tion C{0,j) — Tr{crocr|p} for Hj on a frustrated 11-site 
circle. We assess the quality of the approximation C to 
the exact correlation C by the average relative error 



Error : 



S,|c(o,j)-c(o,j)| 
E,|c(o,j)l 



(7) 



Sliding window is again in very good agreement with the 
exact value for a relatively small window size. For the 
values of Nr < 10 accessible with modest computational 
resources, the replica QBP reproduces the exact corre- 
lation function within a few percents at sufficiently high 
temperatures P G, which is consistent with the system- 
atic error due to the TS decomposition. 

Indeed, both the TS decomposition and the loopy QBP 



4 




FIG. 2: (Color online) Correlations for Hi on a 11-site circle 
at /3 = 6. Exact numerical solution (dash), sliding window 
with i = b (dash-dot), and the replica method for various 
values of Nt (full). Left inset: Error Eq. ([7} vs the Nr for 
different /3. 



contribute to the total error Eq. ([7]) . By brute force com- 
putation, it is possible to determine exactly what frac- 
tion of the error is caused by each of these approxima- 
tions, and in almost all cases we have studied at critical 

both contributions were comparable. Figure [3^ shows 
each contribution to the total error as a function of the 
transverse field g = {g, 0, 0). 

The most successful applications of classical belief 
propagation algorithm are on graphs whose typical loop 
size is very large. This is the case for instance of low den- 
sity parity check codes [1, 0] and spin glasses on Bethe 
lattices Q. Intuitively, one expects a local algorithm 
like belief propagation to be relatively insensitive to the 
large-scale structure of the graph. We expect QBP to 
share this feature, and Figure [8)3 illustrates the effect of 
the loop size on the average relative error of the corre- 
lation function. The oscillatory behavior of the error is 
explained by the frustration present in odd-size circles. 
Save from these oscillations, the results show a global 
improvement as the loop size increases. Errors obtained 
from sliding window (not shown) also show a clear im- 
provement as the loop size increases, but tend to have 
higher errors on even-size loops. 

We have tested QBP on a variety of graphs depicted 
on Fig. |4] a)-d) . The resulting errors in the correlation 
functions are shown in Fig. 3) The computational cost 
is slightly higher for the Heisenberg model because h^y 
do not mutually commute. This restricts the computa- 
tion to lower values of Nr and consequently yields larger 
errors. Modulo this difference, the error is most promi- 
nent for graphs c) and d) which contain loops of size 3. 
In those cases, we found that the QBP algorithm was 
not converging: the magnitude of the errors is consistent 
with the magnitude of the time fluctuations of buv{t) , 
c.f. Eq. ([3]). As expected, the predicted correlation func- 
tion is in much better agreement with its exact value on 




Transverse field g Loop size 



FIG. 3: (Color online) Replica method with Nt — 10. a) 
Different contributions to the error Eq. ([7]) vs transverse field 
for Hi on a 11-sites circle at /3 = 6. b) Error vs loop size, for 
/3 = 1, 2, . . . 10. 
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FIG. 4: (Color online) Error Eq. ^ for Ising and Heisenberg 
models on various loopy graphs using the replica method. For 
the Ising model Nt — 10 and for the Heisenberg model Nt = 
4. [d) is a torus.] 



graphs a) and b) that have only relatively large loops. 



VI. CONCLUSION 

We have numerically characterized the performance of 
the recently proposed QBP algorithm. In the high tem- 
perature phase, both the replica and the sliding window 
QBP algorithms perform remarkably well on a tree with 
modest computational resources, c.f. Fig. [H and offer 
performances similar to TEBD. On loopy graphs, we 
found that the algorithm gives reliable approximations 
when the loop size is large. Most importantly, when the 
results deviated from the exact value, e.g. in the presence 
of small loops, the algorithm did not reach a steady state, 
i.e. the beliefs Eq. ([3]) were highly fluctuating as a func- 
tion of time. This provides an indirect way of assessing 
the validity of the results. 

In , a technique similar to what we have called the 
replica method was used to investigate the phase diagram 
of quantum spin-glasses on Cayley trees. Based on the 
results we have presented, QBP should be suitable to 
study this phase diagram for more general Bethe lattices 
whose typical loop size scales as log \ V\. In the classical 
setting, it has been argued that the physics of random 
Bethe lattices and Cayley trees is greatly different jSj] . We 
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note that the randomness in quenched disordered systems 
should not affect the performences of QBP. In fact, our 
results obtained for random couplings Juv and random 
local fields g are typically in better agreement than the 
ones we have presented. 

Finally, the low temperature phase of these models 
may be accessible using QBP as part of a variational 
approach based on projected entangled-pair states [1^, 
which are a form of 1-bifactor states. QBP can be used 



to approximately compute the reduced state on pairs of 
sites and minimize their energy. We leave the character- 
ization of this approach for a future study. 
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